###
#REPLICATION Script for Golder, Crabtree, and Dhima. 2019.
#Sona Golder, Charles Crabtree, and Kostanca Dhima. 2019. “Legislative Representation and Gender (Bias).” Political Science 71(1): 1-16.
###

# Initial settings
library(foreign)
library(tidyverse)
library(haven)
library(ggridges)
library(scales)
library(estimatr)
library(rio)
library(jtools)
library(ggalt)
library(hrbrthemes)
library(tidyverse)
library(ggplot2)
library(plotMElm)
library(gridExtra)
library(rgdal)
library(sp)
library(interplot)
library(extrafont)
library(ggstance)
library(float)


#set working directory

#read dataset "goldercrabtreedhima_data.dta"
df<- read_dta("goldercrabtreedhima_data.dta")


#set variables as factor variables
df$female_tr <- as.factor(df$female_tr)
df$female_official <- as.factor(df$female_official)
df$national <- as.factor(df$national)
df$day <- as.factor(df$day)


## Responsiveness
# Reply
table(df$reply)
mean(df$reply)*100 # Response rate
mean(df$reply[df$female_tr == 1])*100 # Response rate
mean(df$reply[df$female_tr == 0])*100 # Response rate


# Helpful
table(df$helpful)
mean(df$helpful_bin)*100 # Percent helpful
mean(df$helpful_bin[df$female_tr == 1])*100 # Response rate
mean(df$helpful_bin[df$female_tr == 0])*100 # Response rate


## Estimatation
# Reply
mod <- lm(reply ~ female_tr + female_official + 
            estimated_age + national + day, data = df)
estimatr::commarobust(mod)


# Helpful
mod.help <- lm(helpful_bin ~ female_tr + female_official + 
                 estimated_age + national + day, data = df)
estimatr::commarobust(mod.help)



###
# FIGURE 1
###

# The effect of female_tr does not vary by official gender, legislature level (i.e. professionalization, seniority, etc)
mod.plot <- plot_summs(mod, mod.help, coefs = c("Female" = "female_tr1"),
                       model.names = c("Reply", "Helpful"),
                       color.class = "CUD Bright",
                       ci_level = .95, inner_ci_level = .9, 
                       scale = F, robust = TRUE) +
  labs(title = "") + 
  theme(text = element_text(size = 12),
        axis.text.x = element_text(size = 14),
        axis.text.y = element_text(size = 14),
        axis.title.x = element_text(size = 18),
        axis.title.y = element_text(size = 18),
        plot.title = element_text(size = 22),
        legend.title = element_blank()) +
  labs(x = "Treatment Effect", 
       y = "",
       caption = "Note: Plotted points represent estimated coefficients, thin lines denote 90% confidence intervals, and thick lines indicate 95% confidence intervals.") + 
  ggtitle("Figure 1: Responsiveness",
          subtitle = "Data from 953 NZ elected officials")

mod.plot

ggsave(filename = "../goldercrabtreedhima_replication/Figure1.png", height = 6, width = 10, units = c("in"))



###
# APPENDIX A
###

#Figure 1: Estimated Age of Officials
fill <- "#4271AE"
line <- "#1F3552"

df$age <- df$estimated_age

p.age <- ggplot(df, aes(x = age)) +
  geom_density(fill = fill, colour = line, alpha = 0.6) + 
  scale_x_continuous("Estimated age") +
  scale_y_continuous("Density") +
  ggtitle("Estimated age of officials",
          subtitle = paste0("Data from ", nrow(df), " officials (Mean ", round(mean(df$age, na.rm = T), 3), ", SD: ", round(sd(df$age, na.rm = T), 3), ")")) +
  theme(text = element_text(size = 14),
        axis.text = element_text(size = 14),
        axis.title.x = element_text(size = 16),
        axis.title.y = element_text(size = 16),
        legend.title=element_blank()) +
  theme(plot.title = element_text(size = 17))

p.age

ggsave("../goldercrabtreedhima_replication/appendixa_fig1age.png", width = 8, height = 5, units = "in")


#Figure 2: Official Gender
table(df$gender_official)
df$gender_official <- tools::toTitleCase(df$gender_official)

p.gen <- ggplot(df, aes(gender_official, ..count..)) + 
  geom_bar(aes(fill = gender_official), position = "dodge") + 
  scale_x_discrete("Officials by gender") +
  scale_y_continuous("Count") +
  ggtitle("Officials by gender",
          subtitle = paste0("Data from ", nrow(df), " officials")) +
  theme(text = element_text(size = 14),
        axis.text = element_text(size = 14),
        axis.title.x = element_text(size = 16),
        axis.title.y = element_text(size = 16),
        legend.title = element_blank()) +
  theme(plot.title = element_text(size = 17))

p.gen

ggsave("../goldercrabtreedhima_replication/appendixa_fig2gender.png", width = 8, height = 5, units = "in")



#Figure 3: Officials by level of government
p.nat <- ggplot(df, aes(national, ..count..)) + 
  geom_bar(aes(fill = national), position = "dodge") + 
  scale_x_discrete("Officials by level of government") +
  scale_y_continuous("Count") +
  ggtitle("Officials by level of government",
          subtitle = paste0("Data from ", nrow(df), " officials")) +
  theme(text = element_text(size = 14),
        axis.text = element_text(size = 14),
        axis.title.x = element_text(size = 16),
        axis.title.y = element_text(size = 16),
        legend.title = element_blank()) +
  theme(plot.title = element_text(size = 17))

p.nat

ggsave("../goldercrabtreedhima_replication/appendixa_fig3govtlevel.png", width = 8, height = 5, units = "in")


###
# APPENDIX B
###

#Table 1: LPM Results
#Model 1 - Reply
mod <- lm(reply ~ female_tr + female_official + 
            estimated_age + national + day, data = df)
estimatr::commarobust(mod)

#Model 2 - Helpful
mod.help <- lm(helpful_bin ~ female_tr + female_official + 
                 estimated_age + national + day, data = df)
estimatr::commarobust(mod.help)


###
# APPENDIX C
###

#Table 2: Chi-square Results
#Model 1 - Reply
chisq.test(table(df$reply,
                 df$female_tr))

#Model 2 - Helpful
chisq.test(table(df$helpful_bin,
                 df$female_tr))


#### APPENDIX D
#Table 3: Logit Model Results
#Model 1 - Reply
mod.log <- glm(reply ~ female_tr + female_official + 
                 estimated_age + national, 
               data = df,
               family = binomial())
summary(mod.log)

#Model 2 - Helpful
mod.help.log <- glm(helpful_bin ~ female_tr + female_official + 
                      estimated_age + national, 
                    data = df,
                    family = binomial())
summary(mod.help.log)

###
# APPENDIX E
###

#Figure 2
library(ri2)
declaration <- declare_ra(N = nrow(df), m = as.integer(nrow(df)/2))
df$Y <- df$reply
df$Z <- df$female_tr
ri2_out <- conduct_ri(
  formula = Y ~ Z,
  declaration = declaration,
  sharp_hypothesis = 0,
  data = df
)
summary(ri2_out)
r <- plot(ri2_out)
df$Y <- df$helpful_bin
ri2_out_help <- conduct_ri(
  formula = Y ~ Z,
  declaration = declaration,
  sharp_hypothesis = 0,
  data = df
)
summary(ri2_out_help)
r.help <- plot(ri2_out_help)

grid.arrange(r, r.help, ncol = 2)

###
#END
###


